MIT-CTP 4076 
CERN-PH-TH-2009/175 

Quantum Adiabatic Algorithms, Small Gaps, and Different Paths 

Edward Farhi, 1 Jeffrey Goldstone, 1 David Gosset, 1 
Sam Gutmann, 2 Harvey B. Meyer, 1 ' 3 and Peter Shor 1,4 

1 Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139 
2 Department of Mathematics, Northeastern University, Boston, MA 02115 
3 Physics Department, CERN, 1211 Geneva 23, Switzerland 
^Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139 

O Abstract 

We construct a set of instances of 3SAT which are not solved efficiently using the simplest quantum adi- 
abatic algorithm. These instances are obtained by picking random clauses all consistent with two disparate 
planted solutions and then penalizing one of them with a single additional clause. We argue that by ran- 

i 1 domly modifying the beginning Hamiltonian, one obtains (with substantial probability) an adiabatic path 

that removes this difficulty. This suggests that the quantum adiabatic algorithm should in general be run 
on each instance with many different random paths leading to the problem Hamiltonian. We do not know 
whether this trick will help for a random instance of 3SAT (as opposed to an instance from the particular set 
we consider), especially if the instance has an exponential number of disparate assignments that violate few 
clauses. We use a continuous imaginary time Quantum Monte Carlo algorithm in a novel way to numerically 

(N 

investigate the ground state as well as the first excited state of our system. Our arguments are supplemented 
\q by Quantum Monte Carlo data from simulations with up to 150 spins. 
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I. INTRODUCTION 



Quantum adiabatic algorithms are designed for classical combinatorial optimization problems 
[8]. In the simplest case, such algorithms work by adiabatically evolving in the ground state of 
a system with Hamiltonian H(s) = (1 — s)Hb + sHp that is a function of a parameter s which 
is increased from to 1 as a function of time. Hp is called the beginning Hamiltonian and Hp, 
which is instance dependent, is called the problem Hamiltonian. The minimum (for s £ [0, 1]) 
eigenvalue gap between the ground state and first excited state of H(s) is related to the runtime of 
the adiabatic algorithm. If H(s) has an exponentially small minimum gap then the corresponding 
algorithm is inefficient, whereas a minimum gap which scales inverse polynomially corresponds to 
an efficient quantum adiabatic algorithm. 

Whether or not quantum adiabatic algorithms can be used to solve classically difficult optimiza- 
tion problems efficiently remains to be seen. Some numerical studies have examined the quantum 
adiabatic algorithm on random sets of instances of optimization problems where these sets are 
thought to be difficult for classical algorithms. These studies have reported polynomial scaling of 
the minimum gap out to about 100 bits (6j[l2j[20]. Whether this scaling persists at high bit number 
has recently been called into question |21] . Meanwhile there has been no rigorous analytical result 
that characterizes the performance of the quantum adiabatic algorithms on random instances of 
NP-complete problems. 

Over the years there have been a number of proposed examples which were meant to demonstrate 
failures of the adiabatic algorithm on specific problems. In reference |18J . van Dam et al constructed 
examples intended to foil the adiabatic algorithm, but these examples used a nonlocal cost function. 
Related 3SAT examples of van Dam and Vazirani [19j indeed cannot be solved efficiently using the 
quantum adiabatic algorithm. However it was shown in [5] that such 3SAT instances do not pose a 
problem for the quantum adiabatic algorithm if, having fixed a specific problem Hamiltonian, one 
randomly chooses multiple interpolating paths between the initial and final Hamiltonians and runs 
the adiabatic algorithm once for each random path. Fisher [9] has constructed an interesting but 
specialized example on which the quantum adiabatic algorithm is inefficient. (A later example of 
Reichardt [T7] is based on this.) We do not know if random path change helps here. (It is interesting 
to note that in this case the runtime scales like where n is the number of bits.) The authors 
of reference [22] pointed out that a certain adiabatic algorithm for 3SAT is not efficient. However, 
this was due to a perverse and avoidable nonlocal choice of beginning Hamiltonian Hp [Z] • 

One purpose of this paper is to discuss a different type of challenge to adiabatic optimization 
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Figure 1: Energy levels of the Hamiltonian H'(s) before adding the last term h to the problem Hamiltonian. 
The lower curve coincides with the upper curve only at s = 1. 

which has recently come to light. (See references [2HI] and in a different context |15j.) It seems to 
us that in the history of challenges to adiabatic optimization, this may be the most serious. The 
takeaway message of this work is that a small minimum gap for H (s) can arise when the Hamiltonian 
Hp has features which are seen in the following construction. First suppose that we start with a 
problem Hamiltonian H' p which has two degenerate ground states \z\) and \z2) corresponding to 
bit strings of length n that differ in order n bits. We then expect that the two lowest eigenvalues 
of H'(s) = (1 — s)Hb + sH'p will look something like figure [T] Note that one curve is always below 
the other except at s = 1 even though in the figure they appear to meld because they have the 
same slope at s = 1. Suppose that the upper curve (the first excited state for s near 1) approaches 
the state \z%) and the lower curve approaches the state \z±) as s — > 1. We now form the problem 
Hamiltonian Hp = H' p + h, where h is a term that penalizes the state \z\) but not the state l^)- 
Then, for the two lowest eigenvalues of H(s) = (1 — s)Hp + sHp, we expect to have the situation 
pictured in figure |2j where there is a small gap near s = 1. 

In this paper we discuss a simple way of generating random instances of 3SAT where the cor- 



responding adiabatic Hamiltonian H(s) has the difficulty discussed above. We argue in section III 
that this problem can be overcome (for the set of instances we consider) by randomizing the choice 
of beginning Hamiltonian. 

We present a continuous imaginary time Quantum Monte Carlo algorithm, which is a modifica- 
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Figure 2: Energy levels of the Hamiltoiiian H (s) . There is a tiny gap at s* . 

tion of the Heat Bath algorithm of Krzakala et al |14j , ("Quantum Monte Carlo" is a completely 
classical numerical technique for finding properties of quantum systems and is not a quantum al- 
gorithm.) We use this numerical technique in a novel manner which allows us to investigate both 
the ground state and the first excited state of our Hamiltonian and thereby detect the presence or 
absence of a small gap between them. This differs from the standard application of the Quantum 
Monte Carlo method in that we are able to obtain information about the first excited state using 
a simple procedure which we have validated at low bit number by comparing our results to exact 
numerical diagonalization. These Quantum Monte Carlo simulations support our claim that the 
problem we describe can be overcome using random path change. 

II. PROBLEMATIC INSTANCES 

We now describe the method we use to generate n bit random instances of 3SAT that lead to 
quantum adiabatic Hamiltonians with small minimum gaps. To do this we first generate an instance 
with exactly two satisfying assignments given by the bit strings 111...1 and 000... . Each clause 
c of the 3SAT instance specifies a subset of 3 bits ii(c), 12(c), 13(c) € {l,..,n} and a particular 
assignment W\(c)w2(c)w 3 (c) to those three bits which is disallowed. In order to only generate 
instances which are consistent with the bit strings 111...1 and 000. ..0, we only use clauses for which 

w 1 (c)w 2 (c)w 3 (c) G {100,010,001,110,101,011}. 
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We add such clauses one at a time uniformly at random and stop as soon as these two bit strings 
are the only bit strings which satisfy all of the clauses that have been added. (In practice to check 
whether or not this is the case we use a classical 3SAT solver.) We write m for the total number of 
clauses in the instance. We note that the number of clauses m obtained using this procedure scales 
like n log n. We need this many clauses in order to ensure that each bit is involved in some clause. 
On the other hand, when the number of clauses is 5ralogn the probability of additional satisfying 
assignments goes to zero as n — > oo. 

We now consider the problem Hamiltonian H' p corresponding to this instance, which we define 
to be 



c=l V 

Each term in this sum is 1 if 11(0)12(0)13(0) violates the clause and otherwise. For the beginning 
Hamiltonian we choose 



i=i v y 



(i) 



The two lowest eigenvalues of the Hamiltonian H'(s) = (1 — s)Hp + sH' p will then both approach 
as s — > 1, as in figure [T] The ground state for values of s which are sufficiently close to 1 will 
approach either |000...0) or |lll...l) as s — > 1. Suppose for values of s close to 1 the state that 
approaches |000...0) has lowest energy. Then we add an extra term ho which acts on bits 1, 2 and 
3 and which penalizes this state but not |lll...l) 

and pick 

Hp = Hp + ho . 

(Note that this extra term has a multiplicative factor of |, to avoid any degeneracy of the first 

excited state at s = 1.) In the case where the lowest energy state near s = 1 approaches | 111 1) 

we instead add 

The Hamiltonian H(s) = (1 — s)Hp + sHp is then expected to have a small gap near s = 1 as 
depicted in figure [2] We expect the location s* of the minimum gap to approach s = 1 as n — > 00. 
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Location of the Minimum Gap 

In order to determine the dependence of the location s* of the avoided crossing on the number 
of spins n and the number of clauses m, we consider the perturbative corrections to the energies 
of the states |000...0) and |lll...l) around s = 1. We will show that the low order terms in the 

3 

perturbation series reliably predict a crossing at 1 s* = 1 — 0(^j (^) 3 ) ■ 

With the Hamiltonian constructed in the previous section, one of these states has energy at 
s = 1 (call this the lower state \zi)) and the other state has energy \ at s = 1 (we call this the 
upper state \zjj)). We write our Hamiltonian as 

71 

H(s) = (l- S )- + s 

and we then consider the term — (^p) X^=i if as a perturbation to Hp expanding around s = 1. 

To understand why we trust perturbation theory to predict the location of the near crossing, 
consider a system which is composed of two disconnected sectors, A and B, so the corresponding 
Hamiltonian is of the form 

H A (s) 
H B (s) 

In this situation the generic rule that levels do not cross does not apply and we can easily imagine 
that the two lowest levels look like what we show in figure [3| where the levels actually cross at s*. 

Imagine that low order perturbation theory around s = 1 can be used to get good approximations 
to the ground state energies of Ha(s) and Hb(s) for s near s*, even for s somewhat to the left of s*. 
Then it is possible to accurately predict s*. Our situation is very close to this. We can think of A 
as consisting of the states close to zl in Hamming weight, and B as states close to zjj in Hamming 
weight. Similarly, we view Ha and Hb as the restrictions of H to these sectors. Note that it takes 
n powers of the perturbation to connect |000...0) and |lll...l) and this is why we view A and B as 
essentially disconnected. 

Although figure [3] looks like figure |2j it is figure [2] that depicts the actual situation, where the 
two levels avoid crossing. This true near cross means that the perturbation series in the actual 
theory will diverge very close to s*. However this divergence will only be seen at high order, in fact 
at an order which is proportional to n. The low order terms of the perturbation series in figure [3] 

1 The notation /(n) = Q(g(n)) means that, for n sufficiently large, bg(n) < f(n) < cg(n) for some constants b and 
c. 
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Figure 3: A true energy level crossing can arise from two disconnected sectors. 

are the same as the low order terms of the perturbation series in figure [2] so we can trust low order 
perturbation theory to locate s*. 

(We argue below that as a function of the number of bits, n, s* goes to 1 as n goes to infinity. 
This implies that the radius of convergence of the perturbation theory for the full H(s), expanded 
about s = 1, goes to as n goes to infinity. This fact has no bearing on our argument that low 
order perturbation theory can be used to accurately predict s*.) 

For small values of the parameter = = ^-, the energies of the two states under consideration can be 
expanded as 



El{s) 



Eu(s) 



1 



. n 



. n 



+ 
1 

~ + 



1 



1 



,(2) 



e (4) + 



(3) 
(4) 



It is easy to see that each expansion (inside the square brackets) only contains even powers. Note 

(2) 

that e L is guaranteed to be negative as is always the case for 2nd order perturbation theory of 
the ground state. In addition, since we added a term to penalize the state which had smaller 

(2) (2) 

energy near s = 1 before adding the clause, we expect that e)j < e L . We are interested in the 
behaviour of the difference El(s) —Ejj(s) for randomly generated instances (as generated using the 
prescription of the previous section) as a function of the number of spins n in the limit n — > oo. 
This requires us to further investigate the behaviour of each coefficient e£ anci e[/ as a function 



of n and m. 
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In fact to locate s* we need only to go to second order in perturbation theory. From equation 
[2] we view Hp as the unperturbed Hamiltonian, and — ^ X^I=i °x as ^ ne perturbation, with as 
the expansion parameter. Now <J l x \z) = \z © <§j) so at second order we get 

(2) _ if |(z L eej|4|z L )| 2 



Similarly, 



,(2) 



4 f-* (z L |#p|z L ) - (zjr, © ei|#p|z L © e t 
1 n 1 

4 f-f (z L ©ei|F P |2; L ©ei) 
i=i 



1 n 

z E t 



* i=i 2 (zu ® ei\H P \zu ® £i) 

where the \ in the denominator is (zu\Hp\zu). 

The expected energy penalty incurred when flipping a bit of either zjj or zl is of order ^ since 
each bit is typically involved in O(^) clauses. So the coefficients e£ an d e^f are of order n (^) 
since the energy denominators involved are We now show that their difference is of order 

Write 

V rn ) 



P (2) _ p( 2 ) 



8=1 



where for each i, 



di = - ( 1 h 1 — *— ; r I • (5) 

Recall that f/p = H' p + /i where h is the penalty term from the final clause which acts only on 
the first 3 bits. Therefore, for i = 4,5, n 

1/1 1 

+ 



4\ (zu ® ei\H'p\zu © §i) (z L © ei\H' p \z L © e { 
Our procedure for generating instances is symmetric between the strings 000. ..0 and 111...1 so 
averaging over instances it is clear that the mean of di for i = 4, 5..., n is 0. Thus we expect Ya=a di 
to be (approximately) Gaussian with mean and standard deviation proportional to yJno{d), where 
cr(d) is the standard deviation of each di for i £ {4, 5, n}. To compute u{d) we note that 

!{{zu® ei\H'p\zu © e;) - (z L © ei\H' p \z L © h 



d% 4 V {zL®ei\H' p \z L ®ei){zu®ei\H' P \zu®ei 



Again using the symmetry between all zeros and all ones, we conclude that the numerator is of order 

'C7 



^ and the denominator is of order f^J . Hence we expect a(d) to be 0( (— ) 2 ) . So 
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is of order yfn (5;) 5 ■ We will now locate s* using second order perturbation theory and afterwards 
argue that higher orders do not change the result. Returning to equations [3] and |4j equating the 
two energies at second order we have 



= s* 



1 (l-s* 

■x + 



so 



,( 2 ) _ P ( 2 ) 



l-G(4-(^) 3 ). 



(6) 



The 4th order correction to the energy of the lower state is given by 



,( 4 ) 



( z l\V(—) V\z L ){z L \V—V\z L ) - {zl\V—V— V— V\z L ) 
tip rip tip rip rip 



where 



V 



1 n 

i=l 



1 G ^ ^ CH P (z L © e,i)f Hp{z l © ej) 



and <pL = 1 — Writing rip(z) = (z\Hp\z), this can be expressed as 

16 ^ 

1 E 

1 fi ^ 



16 ^ ^ P (z L © ei)n P (z L © e j )^ P ( 

2 I ® e « ® e j/ 



16 ^ (^p(^l © e;)) 2 n P (z L © ei © %) 
Ajl 1 ^ 1 H P (z L ® ei ® ej] 



Up(zL 



e,: 



Up{z L 



t 16 {n P (z L © ei)) d ^ 16 (H P {z L © ei))" "Hp(z L © ej) {ri P {z L © e< © ej)) 

Now consider the terms in this expression corresponding to indices i,j for which i ^ j and the bits 
i and j do not appear in a clause together. Under these conditions we have 



H P (z L © ii © = %p(z l © gj) + rl P {z L © gj) . 



So we can write 



,( 4 ) 



n 

T- 



r le^p^ffie,)) 3 

1 Up(z l © gj © gj) - ^p(^l © gj) - Up{z L © gj) 

*i cl_tes 16 ( n P( Z L © ^)) 2 «P(^ © © * © g,-)) 



(7) 
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Here the subscript "clausemates" indicates that we sum only over pairs of indices which appear 
together in at least one clause of the 3SAT instance corresponding to Hp. For eff, since the 
unperturbed energy is i we obtain 



e 



n 



(4) 

u fei6(Wp(^ei,)-^ 3 



^ 16 (Hp(z u ®e i )-l) 2 {Hp(.Zu®e j )-l)(Hp(zu®ei®e j )-l) 



«5^j clausemates 



Let's look at the first sum in each of equations [7] and [8] where i goes from 1 to n. Each term as i 
goes from 4 to ra is of order (^) and so the difference of the two sums is e(VnQs) ) ■ The second 
sums (those which are restricted to clausemates) contain of order m terms. Each denominator is 
of order (^) 4 and the numerators are 0(1) since the only contribution to the numerator is from 
clauses in Hp which involve bits i and j together. Separating out the terms where i and/or j are 



1,2 or 3, we conclude that the contribution to — e£* from the clausemate terms is @(^/n [—) 5 ). 
For our instance generation m grows like nlogra so this clausemate contribution is asymptotically 
dominated by the first term which scales like Q(y/n (^) 3 )- So the fourth order contribution to 
the difference of energies Ejj(s) — El(s) is 0(s y / ^(m) 3 (^T^) 4 )■ ^ s * wn i cn we determined 

3 

at second order to be 1 — 0(^yj ^he fourth order contribution to the energy difference is 

O ■ The fourth order corrections can therefore be neglected in determining the location of s*. 

Sixth order and higher contributions to the difference are even smaller. 
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III. FIXING THE PROBLEM BY PATH CHANGE 



In our instance generation we manufactured a small gap by penalizing the planted assignment 
corresponding to the energy eigenstate with the smallest energy near s = 1. Since the slopes of the 
two curves in figure [T] are the same at s = 1, the second derivatives determine which eigenvalue is 

(2) (2) 

smaller near s = 1. After penalization we have e v < e L which is consistent with the near crossing 
in figure [2] Suppose instead that we penalize the assignment corresponding to larger energy in 
figure [TJ Then we expect the situation depicted in figure [4] where no level crossing is induced. 

We imagine that the instances that we manufacture with a small gap as in figure [2] are a model 
for what might be encountered in running the quantum adiabatic algorithm on some instance a 
quantum computer is actually trying to solve. There is a strategy for overcoming this problem. The 
idea is to produce figure [4j with reasonable probability, by randomly modifying the adiabatic path 
which ends at Hp, of course making no use of the properties of the particular instance. For this 
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Figure 4: Energy levels with no avoided crossing near s = 1. Here the second derivative of the upper curve 
B is greater than the second derivative of the lower curve A. 

purpose one could use any random ensemble of paths H(s) such that H(l) = Hp and the ground 
state of H(0) is simple to prepare. However in this paper we only consider randomly changing the 
beginning Hamiltonian. We have made this choice so that we are able to use our Quantum Monte 
Carlo method to numerically verify our arguments. Like other Quantum Monte Carlo methods, 
the method we use does not work when the Hamiltonian has nonzero off-diagonal matrix elements 
with positive sign. 

Starting with an instance where H(s) = (1 — s)Hp + sHp has a tiny gap due to the problem 
discussed above, we now consider a different adiabatic path H(s) = (1 — s)Hp + sHp obtained by 
keeping the same problem Hamiltonian Hp but choosing a different beginning Hamiltonian Hb in 
a random fashion which we will prescribe below. We argue that the small gap near s = 1 is then 
removed with substantial probability, so that by repeating this procedure a constant number of 
times it is possible to find an adiabatic path without a small gap near s = 1. 

The way that we choose a random Hamiltonian Hp is to first draw n random variables a for 
i = 1, 2, n , where each Cj is chosen to be ^ or | with equal probability. We then take 

i=l 

(2) (2) (2) (2) ~ 

We write e)j and e L for the analogous quantities to e L and e\j for the new Hamiltonian H{s) = 
(1 — s)Hb + sHp. The point is that by randomizing Hp in the way we prescribe, there is a 
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substantial probability that one will obtain e\j — e L > , and in that case one expects no avoided 
crossing near s = 1. Write 

n 

e v — e L = c^di 
i=i 

where the {di} are fixed by the instance (and are defined in equation [5J . Since we have fixed the 
problem Hamiltonian Hp, the only random variables appearing in the above equation are the Cj. 
We have cf = | so the mean value of effi — is then 

igr^ = 5 (e «_ ef)<0 . 

But more importantly 

The variance of this difference is 

vJef-~ef) = f^fVar^) 
^ ' i=i 

n 

i=i 

which is 0(n (^;) 3 )- For a fixed instance with a corresponding fixed set {di} the random variable 
c|dj is approximately Gaussian and from its mean and variance we see that the probability 

(2) (2) / \ — 

that e\j — e L is positive and in fact greater than a^fn f^J 2 , for a > 0, is bounded away from 
independent of n. This means that there is a good chance that randomizing Hp turns the situation 
depicted in figure [2] into the situation depicted in figure [1] 

In the case of two planted satisfying assignments with one penalized to produce a small gap 
when the beginning Hamiltonian is Hb of equation [TJ we have shown that a random choice for the 
beginning Hamiltonian Hb of equation [9] can with substantial probability remove the small gap. 
This gives further weight to the idea that when running the quantum adiabatic algorithm on a 
single instance of some optimization problem, the programmer should run the quantum adiabatic 
algorithm repeatedly with different paths ending at Hp 0. 

What if there are k > 2 satisfying assignments and all but one are penalized? 

We have considered instances of 3SAT for which the corresponding problem Hamiltonian has a 
unique ground state and a nondegenerate first excited state that is far from the ground state in 
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Hamming weight. In this case we have shown that the path change strategy succeeds in removing 
a tiny gap. What happens when we have k mutually disparate solutions and we penalize all but 
one of them? We show (under some assumptions) that for large n, path change will succeed after 
a number of tries polynomial in k. We can therefore hope for success if there are k disparate 
assignments which violate few clauses and k scales polynomially with n. On the other hand when k 
is superpolynomial in n we have no reason to be optimistic about the performance of the quantum 
adiabatic algorithm. 

The energies E r (s) for r = 0, 1, .., k — 1 of the k states under consideration can be expanded as 
in equations [3] and [4] Write e^ 1 for the second order corrections in these expansions. If we add 

a small number of clauses to this instance which penalize all the solutions except for a randomly 

(2) (2) 

chosen one which we call zq (forming a new problem Hamiltonian H p ) then the differences e r — 
will not change substantially. So we are interested in the differences 



e ( 2 ) - e (2) 



i=l 



where 



d " <*( {z r e ii\H p \z r e ii) + {z e ei\H P \z e ei)) ' ^ 

and where Zo, Z\, Z2, Zfc— 1 are the k ground states of H p . The adiabatic algorithm applied to 
this instance with problem Hamiltonian H p will succeed if all of the above differences are positive, 
for r = l,...,k — 1. There is a j- chance of this occurring. If they are not all positive, then we 
have encountered an instance which we expect to have a small gap. In this case, we randomize the 
beginning Hamiltonian as in the previous subsection. This produces new values e r for the second 
order energy corrections and the relevant differences are given by 

n 

e~ r ^ — e = ^ ] Cjd r j, (11) 
i=l 

where the random variables q are 5 or | with equal probability (and the d r i are fixed). With this 

fixed problem Hamiltonian, how many times do we have to randomize the beginning Hamiltonian 

(by drawing a random set {q}) to produce algorithmic success? For large n, the random variables 
(2) (2) 

e^. — eg are approximately jointly Gaussian with expectation and variance given by 



s(2) _ s(2) 



5 " 



4 



i=l 



Var (ef > - ef >) = £ (d ri ) 



i=i 



13 



The correlation (the covariance divided by the product of the standard deviations) between any 
two of these random variables is given by 



Corr (ef-ef\~ef)-4 



(2) 



Ej Ki) 2 Ej (drj) 2 



(12) 



We now argue that the above correlation is typically equal to | when n is large. To see this, we can 
expand the above expression using equation [TU] We imagine that the instances have been generated 
so that the law of large numbers applies and that in the limit of large n the quantity 

1 1 



n ^ 



i \ y't 93 &i\Hp\z r © ej) (z, t 
approaches an r independent value and the quantity 

> I — ' 

n 



&i I Hp I %r 



ei) 



1 n 



ei\Hp\z g 



&i | Hp I Z T 



1=1 v H 

approaches a value independent of r and q for r ^ q. Using this in equation 12, we obtain that the 
correlation is equal to ^ (for r ^ q and in the limit of large n) . 

Recall that the problem Hamiltonian has been fixed and therefore the {d r i} are set. We have 

{(2) (2)1 
e T — e > have approximately a joint normal distribution 

with pairwise correlations equal to ^. The probability that a random choice of the {cj} make the 

algorithm succeed is 



Pr 



4 2) -4 2) > , r = l,...,k 



Pr 



Pr 



c?d r i > , r = 1, k — 1 



8=1 



> 



4 Z^i=l "r 



l,...,k 



'T,i=l ( d ri) 2 V E"=l ( d «) 2 

Note that in the last line we have subtracted the mean and divided by the standard deviation. 
Using a standard trick we can rewrite the above probability as (in the limit of large n), 



Pr 



?/;,. 



w > 



.5 V^n , 
4 l^i=l u r 



1. 



,k — l 



(13) 



'E7=i (dn) 2 

where w;o and {ttv} are independent normal random variables with mean zero and variance ^- This 
follows from the fact that the joint normal distribution is specified by its marginal distributions 

(2) 

and pairwise correlations. Essentially w r and wo are the (normalized) changes in the quantities e\. 

(2) 

and eg provided by the randomization. 



At this point, in order to evaluate equation 13 we need to understand the magnitude of the RHS 
of the inequality. Consider typical values for the sum Ei=i dri- As in section [IT] we assume that 
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this sum is Gaussian in the limit of large n with mean and standard deviation approximated by 



(In the argument below, we use some properties of the Gaussian distribution which follow from 
the inequalities 



e~ x dx < I -e~ x dx = — e~ a 
a 2a 

oo ra+1 

e~ x2 dx > / e~ x2 dx > e' {a+1) ' 2 .) 



For the k — 1 approximately Gaussian random variables „ =1 =g with mean and variance 

1, for some constant B it will hold with high probability (in the limit of large n) that 

_5 Y n d ■ 
4 ^ =1 = < f o r r = h •••> * - !• 

E?=i ( d «) 2 

Now the random variable wq satisfies wq < — 2S v / log k with probability which is only polynomially 
small in k. For each w r the probability that w r < —B^logk is polynomially small in k and we use 
the bound 

fe-i 

Pr w r - w > B v / k^k, r = 1, >Pr «) < -2Bylogfcj (l - Pr w r < -B^/\^gk ) . 

r=l 

The first term on the RHS is only polynomially small, and each of the k — 1 terms in the product 
are at least 1 — if B is large enough. So the probability that the adiabatic algorithm succeeds is 
at worst polynomially small in k. By repeating the randomization a polynomial number of times 
we expect to succeed with high probability. 



IV. QUANTUM MONTE CARLO AND NUMERICAL RESULTS 
Continuous Imaginary Time Quantum Monte Carlo 

This section is a review of continous imaginary time Quantum Monte Carlo (which is a classical 
path integral simulation technique for extracting properties of quantum systems |16|). In particular 
we will show how this method can be used to compute thermal expectation values of Hermitian 
operators at inverse temperature (3. We start with a Hamiltonian H which we write as 

H = H + V 

where Ho is diagonal in some known basis {|z)}, and V is purely off diagonal in this basis. We require 
that all the nonzero matrix elements of V are negative. For the Hamiltonian H (s) = (1 — s)Hb+sHp 
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with Hb as in equation [T] we have 



H = sH P + 



(1 — s)n 



V 



Here we include the factor of 1 — s in the definition of V since we are not doing perturbation theory 
in this quantity. 

The partition function can be expanded using the Dyson series as follows 



Tr 



-PH 



Tr 



-PH 



£(-D 



m=0 



rUn rt2 

dt m I dt m —\... / dt l V I (t m )...V I {t l ) 

im=0 Jt m -1=0 JtU=0 



Here we use the notation Vi(t) = e tH °Ve tH °. We can then insert complete sets of states and take 
the trace to obtain the path integral 



Tr 



-PH 



£ E 

m=0 { Zl ,...,z m } 



(-l) m (z 1 \V\z m )(z m \V\z m _ 1 )...{z 2 \V\z 1 ) 



ftm rti „ fl 

dtm / dt m - X ... / A ie -/£*«oW*))* 

tm=0 Jt m -l=0 .7*1=0 



(14) 



In this formula we have used the notation %o{z{t)) = {z(t)\Ho\z(t)) , where the function z(t) is 
defined by 

zi, < t < t\ 

Z 2 , t\ < t < t 2 



z{t) 



zi, t m <t<(3. 



So in particular 

rP 



%o(z(t))dt = (zi|-f/oki)(*i + P -t m ) + (z2\H \z 2 )(h - h) + ...(z m \H \z m )(t m - i m _i) 



t=o 



We view the function z{t) as a path in imaginary time which begins at t = and ends at i = /3. 



Then equation 14 is a sum over paths, where every path is assigned a positive weight according to 
a measure p 



/ 5(P) = (-ir(z 1 |y|z m )(z m |y|z m _ 1 )...(z 2 |y|z 1 )dt 1 ...^ m e-/*=o«oWt))^ 



16 



Note that the fact that p is positive semidefinite follows from our assumption that all matrix 
elements of V are negative. We write 

>-m (15) 

for the normalized distribution over paths. 

We now discuss how one can obtain properties of the quantum system as expectation values with 
respect to the classical measure p. We are interested in ground state properties, but in practice we 
select a large value of f3 and compute thermal expectation values at this inverse temperature. If 
(3 is sufficiently large then these expectation values will agree with the corresponding expectation 
values in the ground state. It is straightforward to show that for any operator A which is diagonal 
in the \z) basis, the thermal expectation value can be written as 

1 -(- I A(z(t))dt) p (16) 



Tr[e~P H ] x /3 j 

where A(z(t)) = (z(t)\A\z(t)) . The notation () p means average with respect to the classical proba- 
bility distribution p over paths. 

This immediately allows us to estimate quantities such as the diagonal part of the Hamiltonian 

Tr\Hne~P H ] 1 ^ 
T r[ e- n 1 = I Ho(z(t))dt) p . (17) 

Another quantity that we find useful in our study is the Hamming weight operator defined by 

i=l V J 



which can also be estimated using equation 16 In order to estimate the thermal expectation value 



of the full Hamiltonian Hq + V we use equation [17] as well as the expression 



Tr[Ve-P H ] m 

Tr[e-f> H ] ~ V P [ ' 



where m is the number of transitions in the path. (Equation 19 is not as simple to derive as equation 



16 ) So by generating paths from the distribution p and then computing averages with respect to 
p, we can evaluate the thermal expectation value of the energy at a given inverse temperature /3 
and by taking (3 sufficiently large we can approximate the ground state energy. 

Generating paths from the distribution p is itself a challenging task. We use a modified version 
of the heat bath algorithm of Krzakala et al |14| which we describe in the appendix. As with other 
Quantum Monte Carlo methods, this algorithm is a Markov Chain Monte Carlo method. In order 
to sample from the distribution p, one defines a Markov Chain over the state space consisting of 
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all paths, where the limiting distribution of the chain is p. To obtain samples from p one starts in 
an initial path Zi n it(t) and then applies some number N equ u of iterations of the Markov Chain. If 
N equ ii is sufficiently large then the distribution of the paths found by further iteration will be close 
to p. We use these subsequent paths to compute averages with respect to p. 

Equilibration of the Quantum Monte Carlo and Identification of Level Crossings 

The discussion in the previous section demonstrates how one can estimate expectation values of 
various quantities in the ground state of a quantum system. The method is "continuous time" so 
there is no discretization error and for fixed /3 quantities can be arbitrarily well approximated with 
enough statistics. We use the Quantum Monte Carlo method in a nonstandard way in order to be 
able to study the lowest two eigenstates of our Hamiltonian (as opposed to just the ground state). 

As described in the previous section, the standard procedure for generating configurations is to 
equilibrate to the distribution p from some initial path Zi n it(t) (we call this the seeded path) by 
applying the Monte Carlo update N equ n times. In order to ensure that one has equilibrated after 
N eq uii Monte Carlo updates, one could for example do simulations with two or more different initial 
paths (seeded paths) and check that the values of Monte Carlo observables appear to converge to 
the same seed independent values. 

For each instance we consider, we run two different Monte Carlo simulations that are seeded with 
two different paths it (t) and zj nit (t). These seeds are paths with no flips in them, corresponding 
to the two states 000. ..0 and 111...1, 



With each seed, we run the Monte Carlo simulation for some total number N tota i of Monte Carlo 
sweeps, taking data every k sweeps. Here a single sweep is defined to be n iterations of our Monte 
Carlo update rule as described in the appendix (where n is the number of spins). We then remove 
the first N equ n Monte Carlo samples, and use the remaining Monte Carlo samples to estimate the 
thermal averages 




000... for all t G [0,0] 



111...1 for all t S [O,0\ . 




Tr[He~P H ] 
Tr[e-P H ] 



(W) 




Tr[e-P H ] 



using equations 17|19 and 18 
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In order to convince the reader that the Monte Carlo algorithm works correctly, and that we 
can obtain information about the first excited state, we show data at 16 bits where exact numerical 
diagonalization is possible. In figures [5] and [6] we show the result of this procedure for a 3SAT 
instance with 16 bits, with problem Hamiltonian H p (before adding the penalty term h, so that the 
levels are degenerate at s = 1 ). This instance has 122 clauses. The inverse temperature is f3 = 150. 
The total number of Monte Carlo sweeps at each value of s is N tota i = 200000, with data taken 
after every fifth sweep (this gave us 40000 data samples). We use N equ u = 2500 samples solely 
for equilibration at each value of s. Note that the Monte Carlo simulations with the two different 
seeds (corresponding to the circles and crosses in the figures) only agree for values of s less than 
roughly 0.4. We interpret this to mean that at these values of s the Quantum Monte Carlo has 
equilibrated to the proper limiting distribution p regardless of the seed. As s increases past 0.4, 
the two simulations abruptly begin to give different results. In this case the simulation seeded with 
zj nit (t) finds a metastable equilibrium of the Markov Chain (i.e not the limiting distribution p), 
which corresponds in this case to the first excited state of the Hamiltonian. For s above 0.4 we can 
see from the comparison with exact diagonalization that the two differently seeded values allow us 
to compute the energies (figure [5| and Hamming weights (figure [6| of the lowest two energy levels 
of our Hamiltonian. 

What is happening here is that for large s, the quantum system can be thought of as consisting 
of two disconnected sectors. One sector consists of states in the z basis with low Hamming weight 
and the other with Hamming weight near n. The Quantum Monte Carlo "equilibrates" in each 
sector depending on the initial seed as can be seen by the smooth data in each sector. (Note that 
data is taken independently at each value of s.) The lower of the cross and circle at each value of s 
in figure [5] is the ground state energy. Of course if the classical algorithm ran for long enough then 
the circles in figure [5] would lie on the crosses for every value of s. 

In figures [7] and [8] we show data taken for the same instance after the penalty clause is added that 
removes the degeneracy at s = 1. In this case the two levels have an avoided crossing as expected 
and which can be seen by exact numerical diagonalization in figure [7] However in the Monte 
Carlo simulation with two seeds there are two essentially disconnected sectors and the two levels 
being tracked do cross. When we see this behaviour in the Monte Carlo simulation, we interpret 
it as compelling evidence that there is a tiny gap in the actual system, which occurs where the 
curves cross. In figure [8] first look at the curves coming from the exact numerical diagonalization. 
The Hamming weight of the ground state decreases and then abruptly rises at the location of the 
minimum gap. The Hamming weight of the first excited state also undergoes an abrupt change 
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at the location of the minimum gap. Not surprisingly the exact diagonalization clearly shows the 
behaviour we expect with the manufactured tiny gap. Now look at the Monte Carlo data which is 
interpreted by first looking at figure [7J In figure [7J the true ground state is always the lower of the 
circles and the crosses. For s below s* ~ 0.423 it is the crosses and for s larger than s* it is the 
circles. Accordingly in figure [8] the Hamming weight of the ground state is tracked by the crosses 
for s below s* and by the circles for s above s*. We can therefore conclude, based only on the 
Monte Carlo data, that the Hamming weight of the ground state changes abruptly. At higher bit 
number we do not have exact numerical diagonalization available but we can still use the Quantum 
Monte Carlo to extract key information. 



Randomizing the Beginning Hamiltonian 

For the 16 spin Hamiltonian depicted above, we now consider randomizing the beginning Hamil- 
tonian Hb as described in section |III| In order to minimize the number of times we run the Monte 
Carlo algorithm (this will be more of an issue at high bit number where simulations are very time 
consuming), we generated many different sets of coefficients {c.;} and calculated the differences 
efj — e L with fixed problem Hamiltonian H' p and beginning Hamiltonians 

H B = 2_^ c 

i=l 



According to our discussion in section III, we expect the avoided crossing near s = 1 to be 



~(2) (2) 

removed for choices of coefficients Cj such that the difference efj — e L > . We made a histogram 

of the eff — e^ 1 (shown in figure [9) and randomly chose three sets of coefficients {cj} such that 
(2) (2) 1 

Sjj — e L J > 2 f° r each. Our analysis predicts that in these cases the crossing will be removed. 
As expected, each of these three sets of coefficients resulted in an adiabatic Hamiltonian with the 
small gap removed, although in one there was another small gap at a smaller value of s. We plot 



the Monte Carlo data for one of these sets in figures 10 and 11 

In section [TT] we argued that we could manufacture a small gap in a 3SAT instance as sketched 
in figure [2] Our 16 bit instance is a concrete example of this and the tiny gap can be seen in figure 
[7J In section |III| we argued that by randomizing the beginning Hamiltonian we should be able to 
produce figure [4] This is what we see concretely at 16 bits in figure 10 We now tell the same story 



at 150 bits using only Monte Carlo data since exact numerical diagonalization is not possible. 
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Data for an Instance of 3SAT with 150 bits 



In addition to validating our method at 16 bits, we studied 3SAT instances with 25, 75 and 150 
bits using our Quantum Monte Carlo simulator. The data from these simulations for the most part 
supported our arguments. In this section we present Monte Carlo data taken for a double plant 
instance of 3SAT with 150 spins and 1783 clauses. In this case the inverse temperature f3 = 300, 
and the total number of Monte Carlo sweeps at each value of s is N total = 100000, with data taken 
every fifth sweep (giving 20000 data samples). The first 2500 data samples at each value of s are 
removed for equilibration. We ran our simulations on a 648 processor SiCortex computer cluster in 
embarrassingly parallel fashion. We used a different processor for each value of s and each value of 
the seed. Data taken at the lower values of s took the longest to accumulate, in some cases more 
than 10 days on a single processor for a single data point. 



Figures [12] and 13 show the energy and the Hamming weight before the penalty clause is added, 



running with two different seeds. In figure 12 for large values of s, the crosses are always below the 



circles and track the ground state which ends at 000... as can be seen in figure 13 We have also 
plotted the second order perturbation theory energies for these two levels, expanding around s = 1. 
Note the good agreement between the second order perturbation theory and the Monte Carlo data 
all the way down to s = 0.35. 

Since the lower curve corresponds to 000... 0, we penalize this assignment by the addition of a 



single clause to form Hp attempting to manufacture a near crossing. In figure 14, the circle data 
is below the cross data for s near 1 but the two curves cross near s = .49. This is seen clearly in 



figure 15 where the energy difference is plotted. From the Monte Carlo data we conclude that H(s) 
has a tiny gap. The location of the avoided crossing is well predicted by second order perturbation 



theory as can be seen in figure 15 



The Monte Carlo data at 150 bits shows that we can make an instance of 3SAT with a tiny 
gap following the procedure outlined in section [IT] We now use the Monte Carlo simulator to show 
that a randomly chosen Hg can alter the Hamiltonian H(s) so that this small gap becomes large. 

(2) (2) 

For the instance at hand we first compute the difference — e L for 100000 randomly chosen 



sets of coefficients. The histogram of these differences is plotted in figure 17 After doing this we 

(2) (2) 1 

randomly selected (from this set) two sets of coefficients such that e\j —e L > k- For both of these 
sets of coefficients we saw that the crossing at s ~ 0.49 was no longer present, although in one case 



there appeared to be a new crossing at a much lower value of s. Figures 18 and 19 show the Monte 
Carlo data for the choice of Hp which does not appear to have any crossing. At 150 bits we see 
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compelling evidence that the story outlined in sections [XT] and 



III 



is true. 



Is s* near 1? 



We argued that s* should go to 1 for large enough n. However our 16 bit example has s* ~ 0.42 

_ (m 

/4 V n , 



and the 150 bit example has s* rj 0.49. Recall from equation 1 1 that s* = 1 — 0(— r- ^ — ^ 4 



At 16 bits with m = 122 we have -Kr (-) 4 = 2.29 and at 150 bits with m = 1783 we have 

3 

(n 1 ) 1 = 1-83- Although asymptotically m is of order nlogn, for these values of n we are not 

3 

yet in the regime where (^) 4 "C 1. 

Even though s* is not near 1 for our instance at 150 bits, we see from figure 15 that second 
order perturbation theory can be used to predict the location of s*. This is because the fourth 
order contribution to the energy difference is quite small. So already at 150 bits we can predict the 
presence or absence of an avoided crossing using second order perturbation theory. 



V. CONCLUSIONS 



We have introduced a new Quantum Monte Carlo technique to analyze the performance of 
quantum adiabatic algorithms for instances of satisfiability. Using seeded configurations (in the 
Monte Carlo simulation) corresponding to disparate low lying states, our technique exposed the 
presence or absence of an exponentially small gap without ever actually computing the gap. 

We used this method to numerically investigate a set of random instances of 3SAT which were 
designed to expose a weakness of the adiabatic algorithm. We confirmed that this weakness can be 
overcome for our set of instances by using path change. Our numerical work and the main part of 
our analysis pertains to instances of 3SAT with 2 planted satisfying assignments. 

We have also considered the scenario where an instance has k satisfying assignments and then 
all but one are penalized with a small number of clauses. However, in the k > 2 case we have 
made certain assumptions which make our analysis possible and we do not know if they apply 
to randomly generated instances. In our scenario the adiabatic algorithm with path change will 
succeed in a number of tries which is polynomially large in k. Here k is fixed and n goes to infinity 
but we take the calculation in section |III| as an indication that the algorithm will succeed when k 
grows polynomially with n. 

Throughout this paper we have assumed that our instances have a unique satisfying assignment. 
In this case we believe that the crucial distinction is whether there are polynomially many or 
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exponentially many low lying mutually disparate states above the unique satisfying assignment. 
If there are polynomially many, we have argued that the adiabatic algorithm with path change 
will succeed, but if there are exponentially many we have no reason to be optimistic about the 
performance of the algorithm. 

When there are exponentially many satisfying assignments our analysis does not apply. This 
situation was considered in recent work by Knysh and Smelyanskiy |13j . 

Our results give further evidence that path change must be considered an integral part of the 
quantum adiabatic algorithm. For any given instance, the algorithm should be run with many 
different randomly selected paths which end at the problem Hamiltonian. As long as the algorithm 
succeeds on at least a polynomially small fraction of the trials, it can be used to solve decision 
problems. 
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Appendix: A modified Version of the Heat Bath algorithm of Krzakala et al 



The authors of |14| give a Quantum Monte Carlo algorithm for spin systems in a transverse field. 
The algorithm we use, which is described in this section, is a modified version of that algorithm. 
The modification which we have made is described in subsection [2] everything else in this section 
constitutes a review of reference [14] . Like other worldline Quantum Monte Carlo techniques, 



the algorithm samples the appropriate probability distribution p (see section IV) over paths in 
imaginary time via Markov Chain Monte Carlo. However this algorithm is only applicable to the 
case where the Hamiltonian is of the form H = Hq + V, where Hq is diagonal in the computational 
basis \z) and 



Ci(7 x 



i=l 



for some set of coefficients {cj} which are all positive. 



For a Hamiltonian of this form, the distribution p over paths (from equation 15 ) is given by 

where in this case a path is specified by an n bit string z\ (call this the starting state) at time t = 
and a sequence of flips which occur in bits labeled ii, ...,i m at times t\, ...,t m (which are ordered), 
where each i r 6 {1, ...,n} and t r £ [0,/3] for r £ {1, ...,m}. Another way of specifying a path is to 
specify the path Pj of each spin j 6 {1, n}. So a path P of the n spin system can be written 
P = (Pi, P2, P n ). For each j £ {1, .., n}, Pj specifies the jth bit of the starting state z\ as well as 
the times at which bit flips occur in bit j . Note that we only need to consider paths which flip each 
bit an even number of times, since only these paths occur with nonzero probability. An example of 
a path for a system with 2 spins is given in figure [20] 
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t=p 




Figure 20: An example of a path P for 2 spins. The red dots indicate bit flips. In this example the starting 
state Z\ = 10 and there are four flips in the first spin and 2 in the second (i.e ii, 12, i±, ie = 1 and z 3 , i 5 = 2). 
The times of these flips are labeled t\ through t§. 

We now define in detail the Markov Chain which has limiting distribution p. It will be useful 
to define pj(Pj\Pi, Pj-i, Pj+i, • ••> Pn) (f° r each j G {l,...,ra}) to be the conditional probability 
distribution of the path of the jth bit, conditioned on the remainder of the path being fixed. 

As in |14j . the update rule for the Monte Carlo algorithm consists of the following 3 steps: 

1. Randomly and uniformly choose a spin j (where j G {1, ...,n}). 

2. Remove all flips in the path which occur in bit j. Also remove the jth bit from the starting 
state z\. This step corresponds to wholly removing the path Pj of spin j. 

3. Draw a new path Pj (starting value and flip times) for spin j from the conditional distribution 

Of course, the nontrivial part of this algorithm is in specifying a procedure which executes step (3) 
in the above. To do this, the authors of reference p3] note that for any index j £ {1, ...,n}, the 
diagonal part of the Hamiltonian can be written as 

H = gj + fj4 (20) 

where gj and fj are operator valued functions of all {cr^} except for a J z . For a given path P, we 
define Fj(z{t)) = (z(t)\fj\z(t)) for t G [0,/3]. The function Fj(z(t)) is then piecewise constant, and 
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for a given path can be written as 

(zi\fj\zi), 0<t<h 

(z2\fj\Z2), h < t < t 2 
\Zm | fj I %m) j tm—1 t <C t m 

^(zi\fj\zi), t m <t<f3. 

Although the above expression involves all the bit flips in the path, the function J : j{z{t)) can 
actually only change value at times where bit flips occur in bits other than the jth bit, since the 
operator fj does not involve a J z . Let us then write 

ho, = i <t< h 
h\, ti < t < t2 



h q , L<t<P = I 



9+1 



In this expression the times t s correspond to times at which bit flips occur in bits other than bit j 
(also, h„ = ho). 

With this notation, the procedure of reference |14| that generates a new path for bit j consists 
of the following: 

1. Compute the value of J-j(z(t)) as a function of imaginary time along the path. 

2. In this step we generate boundary conditions for the path of bit j at times to,ti,...t q . In 
other words we choose q + 1 values sq, s q £ {0, 1} such that at time t r the bit j will be 
set to the value s r in the new path that we are generating. To do this, we sample the values 
so, s g £ {0, 1} for spin j at the times to, t q from their joint distribution, which is given 
by 

(s Q \A q \s q ){s q \A q -x\s q -\)...(s\\A Q \s Q ) 



Z(s ,si, ...,s q \{h , ...,h q },{ti, ...,t q }) 



Tr[A q A q -i...A D 



where 



At = e" 



and Aj = tj+i — U . We also define s q+ i = sq ■ 
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3. Having chosen boundary conditions, we now generate subpaths for bit j on each interval 
[U, U+i] °f length Aj. Such a subpath is specified by a number of flips w and the time offsets 
n, T2, t w G [0, Aj] at which flips occur (note that the starting value of the bit is determined 
by the boundary conditions). The number of flips is restricted to be either even or odd 
depending on the boundary conditions that were chosen for this interval in the previous step. 
In this step, the subpath for each interval is drawn from the distribution 

9i(r, ...,r w ) = l^c.« e --ihiKn-0)-(7 a -r 1 )+( 7S -7 a )-...+(Ai-r„)] d ^ (21) 

(si +1 \Ai\si) 

4. Put all the subpaths together to form a new path Pj for bit j on [0, /?]. 

In section [T] we review the method outlined in p3] for sampling from the distribution 
Z(sq, si, 8 q \ (ho, h q ), (t\, ...,t q )) in step (2) of the above. In section^jwe outline our method for 
sampling from the distribution gi{j\, ...,t w ), which differs from the method suggested in reference 

M- 

1. Generating Boundary Conditions for A Single Spin Path 

The prescription outlined in |14| for generating a set So, s\, s q from the distribution 

Z(s , si, Sq\{h , riq), (ti, ...,t q )) Tr\A A — 1 ^ ] 

is as follows. First generate so £ {0, 1} according to the distribution 

p(*o) = Tr[ ^ 1 .„ i4o] («ol^Vi.--^l-o> • 

(computing these probabilities involves multiplying q two by two matrices). Then generate s± G 
{0, 1} according to 

pOiM = i \ \ \ — 1 a A i v (s \A q A q _ 1 ...A 1 \si)(s 1 \A \s } . 
(so\A q A q -i...A 1 Ao\s ) 

Then generate S2 G {0, 1} from the distribution 

P(S2\S1,S ) = — r(sQ\AqAq- 1 ...A 2 \s2)(s2\A 1 \s 1 ) 

{So\A q A q _i...Ai\Si) 

and so on. Note that this generates the correct distribution since 

P(so)p(si\s )p(s 2 \si, S )...p(s q \Sq-l, Sq) = Z(s , S 1; S q \(h , h q ), (ti, ...,t q )) . 
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2. Algorithm for Sampling from the Single Spin Path Integral with Fixed Boundary Con- 
ditions 



We now present an algorithm which samples from the normalized probability distribution over 
single spin paths S(t) for t G [0, A] (S(t) takes values in {0,1}). Here we only present the case 
of paths with boundary conditions 5(0) = and 5(A) = 1. The other three cases are completely 
analogous. We can parameterize such a path with fixed boundary conditions by the times {n, t w } 
at which S(t) changes value. Note that the number w of such flips is odd due to our choice of 
boundary conditions. The distribution over paths that we aim to sample from is given by 

c tO e -ft[(Tl-0)-(75-Tl) + (7S-7a)-... + (A-T ,)] dri-- _ d7Vo 

It will also be useful for us to make the change of variables from the times (ti,...,t w ) to the waiting 
times (ui, ...,u w ) defined by 

ui = n 

u 3 = T 'j ~ r i-i 3 > 2 

Then the weight assigned to each path is 

cWe -h[u 1 -u2+u 3 -...+\-Y% =1 v,k\du 1 ...du w 



which we can also write as 



c w e -f t i h ( 1 - 2S ( t )) dt du 1 ...du, u 



(l\ e -\[ha z -ca x ]\0} 

The algorithm is as follows: 

1. Start at t=0 in state 5(0) defined by the boundary conditions. Define B\ = 1 — 25(0). Set 
i=l. 

2. Draw the waiting time U{ until the next flip from the distribution 

f( Ui ) = [Vh 2 + c 2 + B i h]e- u ^ VW +^ +B ^ 

If Yl)=i u j > ^ tnen 8° to ste P 3- Otherwise define = —Bi and set i — > i + 1 and repeat 
step 2. 
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3. Take the path you have generated (which will in general be longer than A), and look at the 
segment [0, A]. If this path satisfies the boundary condition at t = A then take this to be the 
generated path. Otherwise, throw away the path and repeat from step (1). 



We now show that this algorithm generates paths from the distribution ( 22 ) . Before conditioning 
on the boundary conditions being satisfied, the probability of generating a sequence of waiting 
times in (ux,Ui +dui), (u2, u 2 + <iu2), ...(u w ,u w -\-du w ) followed by any waiting time li^+i such that 

u w+ i > a - Y^j=i u i is g iven b y 

w 

f(u 1 )f(u 2 )...f(u w )du 1 du 2 ...du w Pmb(u w+ i > A - ^Uj) 
= /(n 1 )/(n 2 )..J(n w )dn 1 dn 2 ...^e-( A -^i n ^[^^ +B -+ l/l ] 

" (f[[^h? + dz + Bth]) e- sr =1 Er= a ^ dui . dUwe -(x-j:j =1 u 3)[ vm^ + sB w+lh] 

Aj =1 / 

c w e -\Vci+h? e - f t i (l-2S(t))h dui _ ^ Uw ^ if w ig eyen 



i+(!) 2 +^(t: 



e 



\Vc 2 +h 2 e - f t=Q (l-2S(t))h dui _ dUw ^ w ig odd 



In the last line we have used the difference of squares formula to simplify consecutive terms: 
[v c 2 + h 2 + h] [v c 2 + h 2 — h] = c 2 . When we condition on fixed boundary conditions (whatever 
they may be), this generates the correct distribution over paths. 
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Figure 5: The discontinuity in the circle data that occurs near s = 0.4 is a Monte Carlo effect that we 
understand. As can be seen from the exact numerical diagonalization there is no true discontinuity in either 
the ground state energy or the first excited state energy. For s greater than 0.4 the Monte Carlo simulation 
is in a metastable equilibrium that corresponds to the first excited state. The true ground state energy at 
each s is always the lower of the circle and the cross at that value of s. 
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Hamming Weight of a 16 Bit Instance with 2 Planted Solutions 



+ Monte Carlo Seeded with 000.. .0 
O Monte Carlo Seeded with 111... 1 

Exact Ground State Hamming Weight 

Exact First Excited State Hamming Weight 




Figure 6: Together with figure [5] we see that the discontinuity in the circles appears in data for both the 
energy and the Hamming weight. This is not indicative of a phase transition in the physical system (as 
evidenced by the smooth curves computed by exact numerical diagonalization), but is purely a result of the 
way in which we use the Monte Carlo method. 
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Energy of a 16 Bit Instance After Adding the Penalty Clause 



+ Monte Carlo Seeded with 000.. .0 
O Monte Carlo Seeded with 111... 1 
Exact Ground State Energy 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

S 



Figure 7: After adding the penalty clause we see that the energy levels have an avoided crossing at s sa 0.423. 
The inset shows exact numerical diagonalization near the avoided crossing where we can resolve a tiny gap. 
The ground state energy is well approximated by the lower of the circle and cross at each value of s. 
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Hamming Weight of a 16 bit Instance After Adding the Penalty Clause 
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Figure 8: In this figure there is a phase transition which occurs near s ss 0.423. We see from the exact 
numerical diagonalization that the Hamming weights of the first excited state and ground state undergo 
abrupt transitions at the point where there is a tiny avoided crossing in figure [7] There is also a jump in 
the Monte Carlo data plotted with circles that occurs before the avoided crossing: this is a Monte Carlo 
effect as discussed earlier and has no physical significance. If we look at the Monte Carlo data in figure 
[7] we conclude that below s « 0.423 the crosses represent the ground state and after this point the circles 
represent the ground state. From the current figure, along with figure [7] we conclude that the Hamming 
weight of the ground state jumps abruptly at s* « 0.423. 
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Histogram of Curvature Differences (16 spins) 




Figure 9: The histogram of e\j' — e L for 1 million different choices of coefficients Cj shows a substantial tail 

(2) (2) 

for which e\j — e L > 0. These sets of coefficients correspond to beginning Hamiltonians Hb for which we 
expect the small gap in figure [7] to be removed. 
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Energy with Choice of HB that Removes a Cross (16 spins) 




+ Monte Carlo Seeded with 000.. .0 
O Monte Carlo Seeded with 111... 1 
- Exact Ground State Energy 
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Figure 10: A random beginning Hamiltonian removes the crossing seen in figure[7] The problem Hamiltonian 
is the same as that in figure [7] The circles are always below (or equal to) the crosses for all values of s. This 
means that the circles track the ground state for all s and we see no sign of a small gap in the Monte Carlo 
data. This is consistent with the displayed exact numerical diagonalization. 
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Hamming Weight with Choice of HB that Removes a Cross (16 spins) 



+ Monte Carlo Seeded with 000.. .0 
O Monte Carlo Seeded with 111... 1 

Exact Ground State Hamming Weight 
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Figure 11: From figure [10] we see that the ground state of the Hamiltonian corresponds to the circles for all 
values of s and the Hamming weight of the circles here goes smoothly to the Hamming weight of the unique 
satisfying assignment. (The jump in the Monte Carlo data corresponding to the crosses is due to the Monte 
Carlo effect discussed earlier.) 
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Energy of a 150 Bit Instance with 2 Planted Solutions 
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Figure 12: The crosses, which represent the Monte Carlo data seeded with 111...1, are below (or equal to) 
the circle data for all values of s. (This is seen more clearly in the inset which shows the positive difference 
between the circle and cross values). We conclude that the crosses track the ground state energy which is 
smoothly varying. The jump in the circle data is a Monte Carlo effect and for s above 0.2 the circles track 
the first excited state. 
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Hamming Weight of a 150 Bit Instance with 2 Planted Solutions 
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Figure 13: Comparing with figure [T2| we see that the Hamming weight for the first excited state and the 
ground state are continous functions of s . We only obtain data for the first excited state for values of s 
larger than s w 0.2 . 
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Figure 14: Adding the penalty clause makes the cross data go above the circle data at s « 0.49. This 
is shown in more detail in figure [15] where we plot the energy difference between the first two levels as a 
function of s. We interpret this to mean that the Hamiltonian H(s) has a tiny gap at s* « 0.49. 
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Energy Difference between the Two Lowest Levels (150 spins) 
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Figure 15: The energy difference, circles minus crosses, from figure[l4]near the value of s where the difference 
is 0. Note that second order perturbation theory does quite well in predicting where the difference goes 
through zero. 
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Hamming Weight of a 150 Bit Instance After Adding the Penalty Clause 
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Figure 16: Looking at figure [14] we see that the ground state is represented by the crosses to the left of 
s s» 0.49 and is represented by the circles after this value of s. Tracking the Hamming weight of the ground 
state, we conclude that it changes abruptly at s w 0.49. 
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Energy with Choice of HB that Removes a Cross (150 spins) 
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Figure 18: A random choice of coefficients such that — ef^ > | gives rise to an H(s) where there is no 
longer an avoided crossing. The circles here correspond to the ground state for all s since the cross data is 
always above (or equal to) the circle data for all s. This can be seen in the inset where we have plotted the 
energy difference, crosses minus circles. The crosses have a Monte Carlo discontinuity near s « 0.2, after 
which they correspond to the first excited state. 
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Hamming Weight with Choice of HB that Removes a Cross (150 spins) 



+ Monte Carlo Seeded with 000.. .0 
O Monte Carlo Seeded with 111... 1 




"++++++- 



" ++HHHHH t H +H-H-H+^H-^^^ ^ 



0.2 



0.3 



0.4 



0.5 



0.6 



l I I l I I l l I I l 'l I I I l l I I I I I l l l I 



0.7 



0.8 



Figure 19: Looking at figure[T8]we see that the ground state corresponds to the circles for all values of s so we 
see here that the Hamming weight of the ground state goes smoothly to its final value as s is increased. We 
take this as further evidence that this choice of Hb would correspond to success for the quantum adiabatic 
algorithm for this instance. 



45 



